#include "parser.h"
#define C_MODE 1
#define G_MODE 0
#define TRAP_CODE 1
#define BE_CODE 2

char *line_buffer;
const char *dels = " \t()\n=,";

int flag_ITER; 
int flag_SPD;
int flag_SPARSE;
int flag_TRAN;
double flag_ITOL;
char *DC_name;
double DC_start,DC_fin,DC_step;
double TRAN_time_step, TRAN_fin_time, TRAN_timestamp;

int G2_id;
int N_id;
int N, m2, size;

void *MNA_G;
cs *mna_g;
void *MNA_C;
cs *mna_c;
size_t cnt = 0;

void *A;
cs* A_sp;
cs *a;

double *b, *x, *RHS;

/*	local functions	*/
void init();
void finalize();
void handle_linear_element(char type);
void handle_command();

void MNA_G_init();
void MNA_G_create();
void MNA_G_destroy();

void conjugate_gradient_solve();
void bi_conjugate_gradient_solve();
void cholesky_solve();
void gauss_solve();

void assign_final_values();
void print_results();

void init()
{
	hashed_elements = g_hash_table_new(g_str_hash, g_str_equal);
	hashed_elements_p = g_hash_table_new(g_str_hash, g_str_equal);
	hashed_elements_m = g_hash_table_new(g_str_hash, g_str_equal);

	hashed_element_values = g_hash_table_new(g_str_hash, g_str_equal);

	hashed_nodes = g_hash_table_new(g_str_hash, g_str_equal);
	hashed_nodes_by_id = g_hash_table_new(g_str_hash, g_str_equal);	
	hashed_nodes_values = g_hash_table_new(g_str_hash, g_str_equal);	
	
	hashed_nodes_g2 = g_hash_table_new(g_str_hash, g_str_equal);
	hashed_nodes_g2_by_id = g_hash_table_new(g_str_hash, g_str_equal);
	hashed_nodes_g2_values = g_hash_table_new(g_str_hash, g_str_equal);

	plot_nodes = g_hash_table_new(g_str_hash, g_str_equal);

	g_hash_table_insert(hashed_nodes, g_strdup("0"), g_strdup("0"));	
	N_id = 1;
	G2_id = 0;

	line_buffer = NULL;
	flag_ITER = 0; 
	flag_SPD = 0;
	flag_SPARSE = 0;
	flag_TRAN = 0;
	flag_ITOL = 0.0;
	DC_step = 0;
	
	TRAN_time_step = -1;
	TRAN_fin_time = -1;
	TRAN_timestamp = 0;

	TRAN_init();
}

void finalize()
{
	free(line_buffer);
}

void mna_fill_dense(gpointer key, gpointer g_value, gpointer user_data)  
{	
	char *name, *plus, *minus, *val, *g_id;

	int p,m,m2_index;
	double value;
	char type;	
	char *lookup_token = (char *)malloc(sizeof(gchar));
	char *elem = (char *)malloc(strlen((gchar*)key));
	strcpy(elem,key);

	int n = N + m2 - 1;
	int mode = user_data;
	
	
	double *MNA;
	if(mode == G_MODE)
		MNA = MNA_G;
	else
		MNA = MNA_C;	
	
	name = strtok(elem,dels);
	plus = strtok(NULL,dels);
	minus = strtok(NULL,dels);
	
	type = name[0];
	
	lookup_token = g_hash_table_lookup(hashed_nodes, g_strdup(plus));
	p = atoi(lookup_token);
	
	lookup_token = g_hash_table_lookup(hashed_nodes, g_strdup(minus));
	m = atoi(lookup_token);
	
	lookup_token = g_hash_table_lookup(hashed_element_values, g_strdup(name));
	value = strtod(lookup_token, NULL);

	
	if(mode == C_MODE && type == 'C'){
		if(p){
			MNA[(p-1)*n + p-1] += value;					
			
			if(m)
				MNA[(p-1)*n + m-1] -= value;
			
			b[p-1] += value;	
		}
		if(m){
			MNA[(m-1)*n + m-1] += value;					
			
			if(p)
				MNA[(m-1)*n + p-1] -= value;
			
			b[m-1] -= value;
		}
	}
	else if(mode == G_MODE)
	{
		if(type == 'R'){
			if(p){
				MNA[(p-1)*n + p-1] += 1.0/value;						
				if(m)
					MNA[(p-1)*n + m-1] -= 1.0/value;		
			}
			if(m){
				MNA[(m-1)*n + m-1] += 1.0/value;					
				if(p)
					MNA[(m-1)*n + p-1] -= 1.0/value;
			}
		}
		else if(type == 'L'){//pigi tasis 0
			type = 'V';
			value = 0.0;
		}
		else if(type == 'I'){
			if(p)
				b[p-1] -= value;
			if(m)
				b[m-1] += value;	
		}
		if(type == 'V'){
		       
			lookup_token = g_hash_table_lookup(hashed_nodes_g2, g_strdup(name));
			m2_index = N + atoi(lookup_token);
	
			if(p){
				MNA[(p-1)*n + m2_index-1] += 1.0;
				MNA[(m2_index-1)*n + p-1] += 1.0;
			}
			if(m){
				MNA[(m-1)*n + m2_index-1] -= 1.0; 
				MNA[(m2_index-1)*n + m-1] -= 1.0;	
			}
			b[m2_index-1] += value;			   	
		}
	}
}

void mna_fill_sparse(gpointer key, gpointer g_value, gpointer user_data) 
{	
	char *name, *plus, *minus, *val, *g_id;

	int p,m,m2_index;
	double value;
	char type;
	char *lookup_token = (char *)malloc(sizeof(gchar));
	char *elem = (char *)malloc(strlen((gchar*)key));
	strcpy(elem,key);

	
	long *Mi, *Mj;
	double *Mx;
	cs *MNA;
	int mode = user_data;
	//printf("mode %d\n",mode);
	
	if(mode == C_MODE)
		MNA = mna_c;
	else
		MNA = mna_g;

	Mi = MNA->i;
	Mj = MNA->p;
	Mx = MNA->x;
	
	name = strtok(elem,dels);
	plus = strtok(NULL,dels);
	minus = strtok(NULL,dels);
	
	type = name[0];
	
	lookup_token = g_hash_table_lookup(hashed_nodes, g_strdup(plus));
	p = atoi(lookup_token);
	
	lookup_token = g_hash_table_lookup(hashed_nodes, g_strdup(minus));
	m = atoi(lookup_token);
	
	lookup_token = g_hash_table_lookup(hashed_element_values, g_strdup(name));
	value = strtod(lookup_token, NULL);

	if(mode == C_MODE && type == 'C'){
		if(p){
			Mi[cnt] = p-1;
			Mj[cnt] = p-1;
			Mx[cnt] = value;
			cnt++;			
			if(m){
				Mi[cnt] = p-1;			
				Mj[cnt] = m-1;
				Mx[cnt] = -value;
				cnt++;			
			}			
			b[p-1] += value;
		}
		if(m){			
			Mi[cnt] = m-1;			
			Mj[cnt] = m-1;
			Mx[cnt] = value;
			cnt++;			
			if(p){
				Mi[cnt] = m-1;			
				Mj[cnt] = p-1;
				Mx[cnt] = -value;
				cnt++;
			}
			b[m-1] -= value;
		}
	}
	else if(mode == G_MODE){
		if(type == 'R'){
			if(p){
				Mi[cnt] = p-1;
				Mj[cnt] = p-1;
				Mx[cnt] = 1.0/value;
				cnt++;			
				if(m){
					Mi[cnt] = p-1;			
					Mj[cnt] = m-1;
					Mx[cnt] = -(1.0/value);
					cnt++;			
				}			
			}
			if(m){			
				Mi[cnt] = m-1;			
				Mj[cnt] = m-1;
				Mx[cnt] = 1.0/value;
				cnt++;			
				if(p){
					Mi[cnt] = m-1;			
					Mj[cnt] = p-1;
					Mx[cnt] = -(1.0/value);
					cnt++;
				}
			}
		}
		else if(type == 'L'){//pigi tasis 0
			type = 'V';
			value = 0.0;
		}
		else if(type == 'I'){
			if(p)
				b[p-1] -= value;
			if(m)
				b[m-1] += value;	
		}
		if(type == 'V'){
		       
			lookup_token = g_hash_table_lookup(hashed_nodes_g2, g_strdup(name));
			m2_index = N + atoi(lookup_token);

			if(p){
				Mi[cnt] = p-1;			
				Mj[cnt] = m2_index-1;
				Mx[cnt] = 1.0;
				cnt++;
	
				Mi[cnt] = m2_index-1;			
				Mj[cnt] = p-1;
				Mx[cnt] = 1.0;
				cnt++;
			}
			if(m){
				Mi[cnt] = m-1;			
				Mj[cnt] = m2_index-1;
				Mx[cnt] = -1.0;
				cnt++;
				Mi[cnt] = m2_index-1;			
				Mj[cnt] = m-1;
				Mx[cnt] = -1.0;
				cnt++;
			}
			b[m2_index-1] += value;
		}
	}
}

void nodes_print(gpointer key, gpointer g_value, gpointer user_data)
{
	printf("(%s)\t", (gchar *)key);
	printf("%s\n", (gchar *)g_value);
}

void to_file(gpointer key, gpointer g_value, gpointer user_data)
{
	char *search;
	search = g_hash_table_lookup(hashed_nodes_values, g_strdup(key));	
	if(search == NULL)
		search = g_hash_table_lookup(hashed_nodes_g2_values, g_strdup(key));	
	char *val = strdup(search);

	FILE *fp = fopen((gchar *)key,"a");
	fprintf(fp, "%gs\t%s\n",TRAN_timestamp, val);

}

void mna_assign(gpointer key, gpointer g_value, gpointer user_data)
{
	char *val = (char *)malloc(sizeof(gchar));	
	val = strdup(key);
	
	int i = atoi(val);
	if((guint)user_data != 1)
		i += N;

	char *final_val = (char *)malloc(sizeof(double));
	sprintf(final_val, "%f", x[i-1]);

	if(user_data == 1)
		g_hash_table_insert(hashed_nodes_values, g_strdup(g_value), g_strdup(final_val));
	else
		g_hash_table_insert(hashed_nodes_g2_values, g_strdup(g_value), g_strdup(final_val));

	free(val);
	free(final_val);
}
void assign_final_values()
{
	g_hash_table_foreach(hashed_nodes_by_id, mna_assign, (gpointer)1);
	g_hash_table_foreach(hashed_nodes_g2_by_id, mna_assign, (gpointer)0);
}

void transient_set_values(gpointer key, gpointer g_value, gpointer user_data) 
{	
	char *name, *plus, *minus, *val, *g_id;

	int k = (guint)user_data;
	int p,m,m2_index;
	double value, new_val;
	char type;
	
	char *lookup_token = (char *)malloc(sizeof(gchar));
	char *elem = (char *)malloc(strlen((gchar*)key));
	strcpy(elem,key);

	name = strtok(elem,dels);
	plus = strtok(NULL,dels);
	minus = strtok(NULL,dels);
	
	type = name[0];
	
	lookup_token = g_hash_table_lookup(hashed_nodes, g_strdup(plus));
	p = atoi(lookup_token);
	
	lookup_token = g_hash_table_lookup(hashed_nodes, g_strdup(minus));
	m = atoi(lookup_token);
	
	lookup_token = g_hash_table_lookup(hashed_element_values, g_strdup(name));
	value = strtod(lookup_token, NULL);
	


	if(type == 'I'){
		if(p){
			TRAN_get_value(name, TRAN_timestamp, &new_val, k);
			b[p-1] = new_val;
		}
		if(m){
			TRAN_get_value(name, TRAN_timestamp, &new_val, k);
			b[m-1] = new_val;
		}		
	}
	else if(type == 'V'){	      
		lookup_token = g_hash_table_lookup(hashed_nodes_g2, g_strdup(name));
		m2_index = N + atoi(lookup_token);
	
		TRAN_get_value(name, TRAN_timestamp, &new_val, k);	
		b[m2_index-1] = new_val;
	}	
}

void print_results()
{
	g_hash_table_foreach(hashed_nodes_values, nodes_print, NULL);
	g_hash_table_foreach(hashed_nodes_g2_values, nodes_print, NULL);
}

void print_to_file()
{
	g_hash_table_foreach(plot_nodes, to_file, NULL);
}

void handle_linear_element(char type)
{	
	/* read from netlist */
	char *data = (char *)malloc(strlen(line_buffer));
	char *trans_data = (char *)malloc(strlen(line_buffer));

	char *name = strtok(line_buffer,dels);
	char *plus = strtok(NULL,dels);
	char *minus = strtok(NULL,dels);
	char *val = strtok(NULL,dels);
	
	/* assign unique IDs, which will be used in the algorithms */
	if((type == 'V') || (type == 'L')){		
		char *g = (char *)g_hash_table_lookup(hashed_nodes_g2,name);
		if(g == NULL){
			g = (char *)malloc(sizeof(int));
			sprintf(g,"%d",G2_id++);
			
			g_hash_table_insert(hashed_nodes_g2,  g_strdup(name), g_strdup(g));
			g_hash_table_insert(hashed_nodes_g2_by_id,  g_strdup(g), g_strdup(name));
		}
	}
		
	char *p = (char *)g_hash_table_lookup(hashed_nodes,plus);
	if(p == NULL){
		p = (char *)malloc(sizeof(int));
		sprintf(p,"%d",N_id++);

		g_hash_table_insert(hashed_nodes,  g_strdup(plus), g_strdup(p));
		g_hash_table_insert(hashed_nodes_by_id,  g_strdup(p), g_strdup(plus));
		g_hash_table_insert(hashed_elements_p,  g_strdup(name), g_strdup(p));
	}

	char *m = (char *)g_hash_table_lookup(hashed_nodes,minus);
	if(m == NULL){
		m = (char *)malloc(sizeof(int));
		sprintf(m,"%d",N_id++);

		g_hash_table_insert(hashed_nodes,  g_strdup(minus), g_strdup(m));
		g_hash_table_insert(hashed_nodes_by_id,  g_strdup(m), g_strdup(minus));
		g_hash_table_insert(hashed_elements_m,  g_strdup(name), g_strdup(m));		
	}

	g_hash_table_insert(hashed_element_values, g_strdup(name), g_strdup(val));	
	sprintf(data,"%s %s %s %s",name, plus, minus, val);
	
	if((type == 'V') || (type == 'I')){

		char *transient_spec = strtok(NULL, dels);
		if(transient_spec != NULL){

			strcpy(trans_data,transient_spec);
			strcat(trans_data," ");
	
			while(1){		
				char *trans_ops = strtok(NULL, dels);

				if(trans_ops != NULL){
					strcat(trans_data,trans_ops);
					strcat(trans_data," ");
				}
				else{
					break;
				}
			}
		}
		else{
			transient_spec = strdup("DC");
			strcpy(trans_data,transient_spec);
			strcat(trans_data," ");
			strcat(trans_data,val);
		}
		TRAN_new_source(name,trans_data);
	}
	/* store for later use(MNA_G) */
	g_hash_table_insert(hashed_elements, g_strdup(data), NULL);	
	free(data);
	free(trans_data);
}

void handle_command()
{
	char *command;
	command = strtok(line_buffer, dels);

	if(!strcmp(command,".OPTIONS")){
		char *option;

		while( (option = strtok(NULL,dels)) != NULL)	
		{	
			if(!strcmp(option,"SPARSE")){
				flag_SPARSE = 1;
			}
			else if(!strcmp(option,"ITER")){
				flag_ITER = 1;
			}
			else if(!strcmp(option,"SPD")){
				flag_SPD = 1;
			}
			else if(!strcmp(option,"ITOL")){
				flag_ITOL = strtod(strtok(NULL, dels),NULL);
			}
			else if(!strcmp(option,"METHOD")){
				command = strtok(NULL, dels);		
				if(!strcmp(command,"BE"))
					flag_TRAN = BE_CODE;
				else
					flag_TRAN = TRAP_CODE;
			}
		}
	}
	else if( !strcmp(command,".PLOT") || !strcmp(command,".PRINT") ){
		char *option;		
		while( (option = strtok(NULL,dels)) != NULL)	//V,ignore it
		{
			option = strtok(NULL, dels); //actual node
			printf("plot %s\n",option);		
			g_hash_table_insert(plot_nodes, g_strdup(option), NULL);
		}
	}		
	else if(!strcmp(command,".DC")){
			
			char *tok;
			tok = strtok(NULL, dels);
			DC_name = strdup(tok);
			
			tok = strtok(NULL,dels);		
			DC_start = strtod(tok,NULL);

			tok = strtok(NULL,dels);		
			DC_fin = strtod(tok,NULL);

			tok = strtok(NULL,dels);		
			DC_step = strtod(tok,NULL);
	
			printf("%s %g %g %g\n", DC_name, DC_start, DC_fin, DC_step);
	}

	else if(!strcmp(command,".TRAN")){
		TRAN_time_step = strtod(strtok(NULL, dels),NULL);		
		TRAN_fin_time = strtod(strtok(NULL, dels),NULL);
	}
}

void MNA_G_init()
{
	N = (int)g_hash_table_size(hashed_nodes);
	m2 = (int)g_hash_table_size(hashed_nodes_g2);
	size = N - 1 + m2;

	if(flag_SPARSE){
		printf("sparse ");		
		(cs *)MNA_G;
		(cs *)MNA_C;
		(cs *)A;
	}
	else{
		printf("dense ");		
		(double *)MNA_G;
		(double *)MNA_C;
		(double *)A;
	}
}

void MNA_G_create()
{
	if(flag_SPARSE){
		mna_create_sparse(&mna_g, &b, N, m2);
		g_hash_table_foreach(hashed_elements, mna_fill_sparse, G_MODE);
		MNA_G = cs_compress(mna_g);
		cs_dupl(MNA_G);
		cs_dropzeros(MNA_G);
	}
	else{
		mna_create_dense(&MNA_G, &b, N, m2);
		g_hash_table_foreach(hashed_elements, mna_fill_dense, G_MODE);
	}
}

void MNA_C_create()
{
	cnt = 0;
	if(flag_SPARSE){
		mna_create_sparse(&mna_c, &b, N, m2);	
		g_hash_table_foreach(hashed_elements, mna_fill_sparse, C_MODE);		
		MNA_C = cs_compress(mna_c);		
		cs_dupl(MNA_C);
		cs_dropzeros(MNA_C);
		//cs_print(MNA_C,1);
	}
	else{
		mna_create_dense(&MNA_C, &b, N, m2);
		g_hash_table_foreach(hashed_elements, mna_fill_dense, C_MODE);
	}
}

void MNA_destroy()
{
	if(flag_SPARSE){
		cs_spfree(MNA_G);
		cs_spfree(MNA_C);
	}
	else{
		free(MNA_G);
		free(MNA_C);
	}
}

void gauss_solve()
{

	printf("Gauss LU\n");	
	
	if(flag_SPARSE)
		ca_gauss_sparse(A, RHS, size, x);
	else
		ca_gauss_dense(A, RHS, x, size);
}

void cholesky_solve()
{
	printf("Cholesky\n");
	if(flag_SPARSE)
		ca_chol_sparse(A, RHS, size, x);
	else
		ca_chol_dense(A, RHS, size, x);
}

void conjugate_gradient_solve()
{
	printf("CG\n");
	if(flag_SPARSE)
		ca_conjugate_gradient_sparse(A, RHS, size, x, flag_ITOL);
	else
		ca_conjugate_gradient_dense(A, RHS, size, x, flag_ITOL);
}

void bi_conjugate_gradient_solve()
{
	printf("BiCG\n");

	if(flag_SPARSE)
		ca_bi_conjugate_gradient_sparse(A, RHS, size, x, flag_ITOL);
	else
		ca_bi_conjugate_gradient_dense(A, RHS, size, x, flag_ITOL);
}

void normal_mode()
{
	if (flag_SPD && flag_ITER)
		conjugate_gradient_solve();
	else if (flag_ITER)
		bi_conjugate_gradient_solve();
	else if (flag_SPD)
		cholesky_solve();
	else
		gauss_solve();
}

void ca_parse(char *netlist)
{
	if (netlist == NULL){
		printf("Error: You haven't entered a netlist!\n");
		exit(42);
	}
	/*	Initializations	*/	
	init();
	
	FILE *input_stream;	
	if ((input_stream = fopen( netlist, "r")) == NULL){
		printf("Error: no such file (%s)\n", netlist);
		finalize();
		exit(42);
	}
	
	/*	Read netlist line by line	*/	
	size_t len = 0;		ssize_t read = 0;	char type;
	
	while ((read = getline(&line_buffer, &len, input_stream)) != -1)
	{
		type = toupper(line_buffer[0]);
		
		if((type == 'V')||(type == 'L')||(type == 'I')||(type == 'R')||(type == 'C'))
			handle_linear_element(type);
		else if (type == '.')
			handle_command();
	}

	MNA_G_init();
	MNA_G_create();

	double e[size];
	double e_prev[size];	
	int i;	
	for(i = 0; i < size; ++i)
		e_prev[i] = e[i] = b[i];

	int m,k;

	A = MNA_G;
	RHS = b;
	x = (double *)malloc(size*sizeof(double));	
	
	
	
	//DC
	printf("DC:\n");
	normal_mode();
	assign_final_values();
	print_results();
	
	
	if(TRAN_time_step != -1){	//transient mode
		
		printf("TRAN:\n");
		
		MNA_C_create();	
	
		//getchar();
		RHS =  malloc(size*sizeof(double));		
		if(flag_SPARSE){
			cs_spfree(A_sp);
			A_sp = cs_spalloc (size, size, size*size, 1, 1);
			A = cs_compress(A_sp);
			cs_dupl(A);
			cs_dropzeros(A);
		}
		else{		
			A = malloc(size*size*sizeof(double));
		}
		
		//cs_print(A,1);
		m = floor(TRAN_fin_time/TRAN_time_step);
		for(TRAN_timestamp = 0, k = 0; k < m; k++, TRAN_timestamp += TRAN_time_step){	
			printf("time:%g\n",TRAN_timestamp);
			
			for(i = 0; i < size; ++i){				
				e_prev[i] = e[i];
				e[i] = b[i];
			}
			
			//new round resets			
			memset(b,0,size*sizeof(double));
			memset(RHS,0,size*sizeof(double));
			if(!flag_SPARSE)
				memset(A,0,size*size*sizeof(double));

			//fill new b
			g_hash_table_foreach(hashed_elements, transient_set_values, (gint)k);			
					
			if(flag_TRAN == BE_CODE){
				printf("BE\n");		
				if(flag_SPARSE)
					TRAN_backward_euler_sparse(MNA_G, MNA_C, b, x, size, TRAN_time_step, &A, RHS);
				else
					TRAN_backward_euler(MNA_G, MNA_C, b, x, size, TRAN_time_step, A, RHS);
			}
			else{
				printf("TRAP\n");
				if(flag_SPARSE)	
					TRAN_trapezoidial_sparse(MNA_G, MNA_C, e, e_prev, x, size, TRAN_time_step, &A, RHS);
				else
					TRAN_trapezoidial(MNA_G, MNA_C, e, e_prev, x, size, TRAN_time_step, A, RHS);
			}
			//cs_print(A,1);

			memset(x,0,size*sizeof(double));
			
			normal_mode();
			assign_final_values();
			print_results();
			//print_to_file();
			getchar();
		}
		
	}
	else if(DC_step){	//DC sweep mode
		char *lookup;			
		int p,m,m2_index;	
		char type = DC_name[0];

		lookup = g_hash_table_lookup(hashed_elements_p,  g_strdup(DC_name));
		if(lookup == NULL){
			lookup = g_hash_table_lookup(hashed_nodes_g2,  g_strdup(DC_name));
			m2_index = atoi(lookup) + N;
		}
		else{
			p = atoi(lookup);
			lookup = g_hash_table_lookup(hashed_elements_m,  g_strdup(DC_name));
			m = atoi(lookup);
		}
		double val;
		for(val = DC_start; val <= DC_fin; val += DC_step){	
			printf("val:%g\n",val);

			//fill new b
			if(type == 'I'){
				if(p)
					b[p-1] = val;
				if(m)
					b[m-1] = val;	
			}
			else if(type == 'V'){
				b[m2_index-1] = val;	
			}
			memset(x,0,size*sizeof(double));
			
			normal_mode();
			assign_final_values();
			print_results();
			getchar();
		}
	}
	
	printf("%g\n",TRAN_timestamp);
	print_results();
	print_to_file();
	
	/*Finalize */	
	MNA_destroy();
	finalize();
	fclose(input_stream);
}
